Strongly coupled plasma with electric and magnetic charges 
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A number of theoretical and lattice results lead us to believe that Quark-Gluon Plasma not too far 
from Tc contains not only electrically charged quasiparticles - quarks and gluons - but magnetically 
charged ones - monopoles and dyons - as well. Although binary systems like charge-monopole 
and charge-dyon were considered in details before in both classical and quantum settings, it is 
the first study of coexisting electric and magnetic particles in many-body context. We perform 
Molecular Dynamics study of strongly coupled plasmas with ~ 1000 particles and different fraction 
of magnetic charges. Correlation functions and Kubo formulae lead to such transport properties as 
diffusion constant, shear viscosity and electric conductivity: we compare the first two with empirical 
data from RHIC experiments as well as results from AdS/CFT correspondence. We also study a 
number of collective excitations in these systems. 



I. INTRODUCTION 

This paper has two main goals. One is to introduce 
a new view of finite T — fi QCD based on a competition 
between electrically and magnetically charged quasipar- 
ticles (to be referred to as EQPs and MQPs below). The 
second is to begin quantitative studies of many-body ef- 
fects based on these ideas, starting with the simplest ap- 
proach possible, namely by use of classical mechanics and 
basic forces acting between them. 

The motivations and some details related to those 
quasiparticles will be explained in the rest of the Intro- 
duction: but before doing so let us start with a short 
summary of the picture proposed. It is different from 
the traditional approach, which puts confinement phe- 
nomenon at the center of the discussion, dividing the 
temperature regimes into two basic phases: (i) confined 
or hadronic phase at T < Tc, and (ii) deconfined or 
quark-gluon plasma (QGP) phase at T > Tc- 

We, on the other hand, focus on the competition of 
EQPs and MQPs and divide the phase diagram differ- 
ently, into (i) the "magnetically dominated" region at 
T < Te=m and (ii) "electrically dominated" one at 
T > Te=m- In our opinion, the key aspect of the physics 
involved is the coupling strength of both interactions. So, 
a divider is some E-M equilibrium region at intermedi- 
ate T-/i. Since it does not correspond to a singular line, 
one can define it in various ways"'^: the most direct one 
is to use a condition that electric (e) and magnetic (g) 
couplings are equal 



e^/hc = g^jhc = 1/2 



(1) 



The last equality follows from the celebrated Dirac quan- 
tization condition [1] 



with n being an integer, put to 1 from now on. 

Besides equal couplings, the equilibrium region is also 
presumably characterized by comparable densities as well 
as masses of both electric and magnetic quasiparticles^. 

The "magnetic-dominated" low-T (and low-/.t) region 
(i) can in turn be subdivided into the confining part (i- 
a) in which electric field is confined into quantized flux 
tubes surrounded by the condensate of MQPs, forming 
t'Hooft-Mandelstamm "dual superconductor" [2], and a 
new "postconfinement" region (i-b) at < T < Te=m 
in which EQPs are still strongly coupled (correlated) and 
still connected by the electric flux tubes. We believe 
this picture better corresponds to a situation in which 
string-related physics is by no means terminated at T = 
Tc'. rather it is at its maximum there. Then if leaving 
this "magnetic-dominated" region and passing through 
the equilibrium region by increase of T and/or /i, we 
enter either the high-T "electric-dominated" QGP or a 
color-electric superconductor at high-/i replacing the dual 
superconductor (with diquark condensate taking place of 
monopole condensate). 

A phase diagram explaining this viewpoint pictorially 
is shown in Fig.l. 

The paper is structured as follows. We start with a 
mini-review of the subject, including RHIC phenomenol- 
ogy and lattice data relevant for this picture. Then we 
start discussing classical dynamics of electric and mag- 
netic charges interacting with each other. We briefly 
summarize what is known about the two-body problems 
in section II, and proceed to (idealized) 3-body problem, 
namely a motion of a magnetically charged object in a 
field of a static electric dipole, and discuss fiux tubes 
in classical plasmas. We then proceed to Molecular Dy- 
namics (MD) simulations. The main parameters of the 
model are (i) the ratio of MQPs/EQPs concentration and 



eg n 
Yc~ 2 



(2) 



^Another possibility discussed in section I A is the curve of 
marginal stability for a gluon. 



■^Let us however remind the reader that the E-M duality is 
of course not exact, in particular EQPs are gluons and quarks 
with spin 1 and 1/2 while MQPs are spherically symmetric 
"hedgehogs" without any spin. 
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(ii) the ratio of magnetic-to-electric coupling g/e. The 
main issues discussed are how the transport properties 
(in particular the shear viscosity) of the plasma depend 
on them. More specifically, the issue is whether admix- 
ture of weaker-coupled MQPs increases or decreases it. 
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W2(s,pt,M„y,6, A) =< cos{2(t)) > 



(3) 
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FIG. 1. (color online) A schematic phase diagram on a 
( "compactified" ) plane of temperature and baryonic chemical 
potential T — fj,. The (blue) shaded region shows "magneti- 
cally dominated" region g < e, which includes the e-confined 
hadronic phase as well as "postconfined" part of the QGP 
domain. Light region includes "electrically dominated" part 
of QGP and also color superconductivity (CS) region, which 
has e-charged diquark condensates and therefore obviously 
m-confined. The dashed line called "e=g line" is the line of 
electric-magnetic equilibrium. The solid lines indicate true 
phase transitions, while the dash-dotted line is a deconfine- 
ment cross-over line. 



A. Strongly coupled Quark-Gluon plasma in heavy 
ion collisions 

A realization [3,4] that QGP at RHIC is not a weakly 
coupled gas but rather a strongly coupled liquid has led 
to a paradigm shift in the field. It was extensively de- 
bated at the "discovery" BNL workshop in 2004 [5] (at 
which the abbreviation sQGP was established) and mul- 
tiple other meetings since. 

Collective flows, related with explosive behavior of hot 
matter, were observed at RHIC and studied in detail: the 
conclusion is that they are reproduced by the ideal hydro- 
dynamics remarkably well. Indeed, although these flows 
affect different secondaries differently, yet their spectra 
are in quantitative agreement with the data for all of 
them, from tt to fl~ . At non-zero impact parameter 
the original excited system is deformed in the transverse 
plane, creating the so called elliptic flow described by 



where (/> is the azimuthal angle and the others stand 
for the collision energy, transverse momentum, particle 
mass, rapidity, centrality and system size. Hydrodynam- 
ics explains all of those dependence, for about 99% of the 
particles'^. 

Naturally, theorists want to understand the nature of 
this behavior by looking at other fields of physics which 
have prior experiences with liquid-like plasmas. One of 
them is related with the so called AdS/CFT correspon- 
dence between strongly coupled A/'=4 supersymmetric 
Yang-Mills theory (a relative of QCD) to weakly coupled 
string theory in Anti-de-Sitter space (AdS) in classical 
SUGRA regime. We will not discuss it in this work: for 
a recent brief summary of the results and references see 
e.g. [6]. 

Zahed and one of us [4] argued that marginally bound 
states create resonances which can strongly enhance 
transport cross section. Similar phenomenon does hap- 
pen for ultracold trapped atoms, due to Feshbach-type 
resonances at which the binary scattering length a oo, 
which was indeed shown to lead to a near-perfect liquid, 
van Hees, Greco and Rapp [7] studied qc resonances, and 
found enhancement of charm stopping. 

Combining lattice data on quasiparticle masses and in- 
terparticle potentials, one finds a lot of quark and gluon 
bound states [8,9] which contribute to thermodynami- 
cal quantities and help explain the "pressure puzzle" [8] , 
an apparent contradiction between heavy quasiparticles 
near Tc and rather large pressure. The magnetic sec- 
tor discussed in this paper provides another contribution, 
that of MQPs (monopoles and dyons) , which will help to 
resolve the pressure puzzle. 

A very interesting issue is related with counting^ of the 
bound states of all quasiparticles. Here the central notion 
is that of curves of marginal stability (CMS), which are 
not thermodynamic singularities but lines indicating a 
significant change of physics where a switch from one 
language to another (like E^M ) is appropriate or even 
mandatory. 

Let us mention one example related with quite interest- 
ing "metamorphosis" discussed in literature, in the con- 
text of 7V=2 SUSY theories. The CMS in question is 
related with the following reaction 



gluon monopole + dyon 



(4) 



in which the r.h.s. system is magnetically bound pair 
(obviously with zero total magnetic charge). The curve 
itself is defined by the equality of thresholds. 



M(gluon) — M{dyon) + M (monopole) 



(5) 



^The remaining ^ 1% resigning at larger transverse mo- 
menta pt > 2GeV are influenced by hard processes and jets. 
*And prevention of the double counting. 
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As discussed in details by Ritz ct al [10] , inside the re- 
gion surrounded by CMS (in which the state in the r.h.s. 
is lighter than the gluon) even a notion of a gluon as 
a separate state does not exist, and using the "magnetic 
language" (r.h.s.) for its description becomes mandatory. 



B. Classical Molecular Dynamics for non-Abelian 
plasmas 

Another direction, pioneered by Gelman, Zahed and 
one of us [11], is to use experience of classical strongly 
coupled electromagnetic plasma. Their model for the de- 
scription of strongly interacting quark and gluon quasi- 
particles as a classical and non-rclativistic Non-Abelian 
Coulomb gas. The sign and strength of the inter-particle 
interactions are fixed by the scalar product of their classi- 
cal color vectors subject to Wong's equations. The EoM 
for the phase space coordinates follow from the usual 
Poisson brackets: 

= S"^-S^06,, {Ql„ 0^,} ^ r"-^ Ql, (6) 

For the color coordinates they are classical analogue of 
the SU(Nc) color commutators, with /"^"^ the structure 
constants of the color group. The classical color vec- 
tors are all adjoint vectors with a = 1...{N^ — 1). For 
the non-Abelian group SU(2) those are 3d vectors on a 
unit sphere, for SU(3) there are 8 dimensions minus 2 
Casimirs=6 d.o.f.^. 

This cQGP model was studied using Molecular Dy- 
namics (MD), the equations of motion were solved nu- 
merically for n ~ 100 particles. It also displays a 
number of phases as the Coulomb coupling is increased 
ranging from a gas, to a liquid, to a crystal with anti- 
ferromagnetic-like color ordering. There is no place for 
details here, let us only mention that important transport 
properties like diffusion and viscosity vs coupling, note 
how different and nontrivial they are. When extrapo- 
lated to the sQGP suggest that the phase is liquid-like, 
with a diffusion constant D « G.l/T and a bulk viscosity 
to entropy density ratio r}/sK. 1/3. The second paper of 
the same group [11] discussed the energy and the screen- 
ing at F > 1, finding large deviations from the Debye 
theory. 

The first study combining classical MD with quantum 
treatment of the color degrees of freedom has been at- 
tempted by the Budapest group [12]. 



^Although color EoM does not look like the usual canonical 
relations between coordinates and momenta, they actually are 
pairs of conjugated variables, as can be shown via the so called 
Darboux parametrization, see [11] for details. 



C. Electric-magnetic dualities in supersymmetric 
theories 

Progress in supersymmctryc (SUSY) Quantum Field 
Theories was originally stimulated by a desire to get rid of 
perturbative divergencies and solve the so called hierar- 
chy problems. However in the last 2 decades it went much 
further than just guesses of possible dynamics at super- 
high energies. A fascinating array of nonperturbative 
phenomena have been discovered in this context, making 
them into an excellent theoretical laboratory. However 
we think that their relevance to QCD-like theories are 
neither understood not explored in a sufficient depth yet. 

Studies of instantons in these theories have resulted 
in exact beta functions [13] and other tools, which have 
allowed Seiberg to get quite complete picture of the phase 
structure of A/' =1 SUSY gauge theories [14]. 

This was enhanced in the context of M =2 SUSY gauge 
theories by Seiberg and Witten [15], who were able to 
show how physical content of the theory changes as a 
function of Higgs VeVs (in a "moduli space" of possi- 
ble vacua). Singularities in moduli space were identified 
with the phase transitions, in which one of the MQPs gets 
massless. Seiberg and Witten have found a fascinating 
set of dualities, explaining where and how a transition 
from one language to another (e.g. from "electric" to 
"magnetic" to "dyonic" ones) can explain what is hap- 
pening at the corresponding part of the moduli space, in 
the simplest and most natural way. 

One lesson from those works, which is most important 
for us, is what happens with the strength of electric e and 
magnetic coupling g near the phase transition. As e.g. 
monopoles gets light and even massless at some point, 
the "Landau zero charge" in the IR is enforced by the 
[/(I) beta function of the magnetic QEDs, making them 
weakly coupled in IR, g « 1. The Dirac quantization 
(2) therefore demands that the electric coupling must get 
large e >> 1, enforcing the "strongly coupled" electric 
sector in this region. 

Since two pillars of this argument - U{V) beta function 
and Dirac quantization - do not depend on supcrsymme- 
try or any other details of the SW theory, we therefore 
now propose it to be a generic phenomenon. We thus con- 
jecture it to be also true near the QCD deconfinement 
transition T w Tc, explaining why phenomenologically 
we see a strong coupling regime there. 

The high-T limit, on the other hand, is similar to large- 
VEV domain of moduli space: here the SU{N) asymp- 
totic freedom in UV plus screening makes the electric 
charge small. Thus here MQPs are heavy and strongly 
coupled. 

D. Lessons from lattice gauge theory 

Static potentials. One of the principal reasons we 
proposed to change the traditional viewpoint of putting 
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confinement at the center, can be explained using lattice 
data on the T— dependence of the so called "static poten- 
tials" . The traditional reasoning points to the free energy 
F{r, T) associated with static quark pair separated by a 
distance r, and defines the deconfinement as the disap- 
pearance of a (linearly) growing "string" term in it, so 
that at T > Tc there is a finite limit of the free energy 
at large distances, F{oo,T). This phenomenon has often 
been referred to as a "melting of the confining string" at 
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3 T/T, 



(as TSoc{T)) de- 
rived from the free energy of two static quarks separated by 
large distance, in 2-flavor QCD according to [17]. 

However, as explained by Polyakov nearly 3 decades 
ago [16], the string actually should not disappear at Te- 
at this point its energy gets instead compensated by the 
entropy term so that the free energy F ^ U — TS van- 
ishes. As detailed lattice studies revealed, in fact the en- 
ergy and entropy associated with a static quark pair are 
strongly peaked exactly at T « Tc, see Fig. 2. The poten- 
tial energy is really huge there, reaching about 4 GeV(!), 
while the associated entropy reaches equally impressive 
value of about 20. Nothing like that can be explained 
on the basis of Debye-screened weakly coupled gas of 
EQPs - the usual picture of QGP until few years ago. 
We think that the explanation of such large energy and 
huge number ~ exp{20) of occupied states can only be 
obtained if several correlated quasiparticles are bound to 
heavy charges, presumably in the form of gluonic chains 
or "polymers" [9] conducting the electric flux from one 



charge to another. 

Therefore, the "deconfinement" seen in disappearing 
linear term in free energy is actually restricted to static 
(or adiabatically slowly moving) charges, while for finite- 
frequency motion of light or even heavy (charmed) quarks 
one still should find mesonic bound states even in the 
deconfined phase [8] . Lattice studied of light quark and 
charmonium states [18] found that they indeed persist till 
T w 2Tc: this conclusion was dramatically confirmed by 
experimental discovery that J/ip suppression at RHIC is 
smaller than expected and is consistent with a new view, 
that J/V' is not melting at RHIC (where T < 2Tc). 

One set of well-known lattice studies have tried to an- 
swer the following questions: Is a "dual superconductor" 
picture consistent with what is observed on the lattice? 
In particular, is the shape and field distribution inside the 
confining strings in agreement with that in the Abrikosov 
flux tube of a superconductor (Abelian Higgs model)? As 
one can read e.g. in [19], the answer seems to be a def- 
inite yes. Can one define in some way monopoles and 
their paths, and are those (in average) consistent with 
dual Maxwell equations? As one can read in e.g. [20], 
the answer seems to be also yes. 

Unfortunately, those studies (as reviewed in [21]) were 
mostly concentrated in the vacuum T = 0, while we are 
interested by the deconfined plasma T > Tc- Is there any 
general reason to think that MQPs play an important 
role here as well? The most important argument^ is the 
persistence of static magnetic screening at all T, up to 
infinitely high T. 

Screening. Although static magnetic screening was 
shown to be absent in perturbative diagrams [22], it has 
been conjectured by Polyakov [16] to appear nonpertur- 
batively at the "magnetic scale" which at high T is 



A 



M 



e^T 



(7) 



The magnetic screening mass and monopole density 
should thus be 



Mm = Cm A 



3 

M 



(8) 



with some numerical constants Cm, Cn- 

To illustrate current lattice results, we show the T- 
dependence of the electric and magnetic screening masses 
calculated by Nakamura et al [23], see Fig. 3. Note that 
electric mass is larger than magnetic one at high T , but 
vanishes at Tc (because here electric objects gets too 
heavy and effectively disappear). The magnetic screen- 
ing mass however grows toward Tc, which is consistent 
with its scaling estimate 



^Note a principle difference with all electromagnetic plas- 
mas, which have no magnetic screening at all. For example, 
solar magnetic fiux tubes are extended for a millions of km, 
with unimpeded flux. 
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Ml 



{e'Tf 



(9) 



(Another estimate of the magnetic screening can be done 
in the dual language as 



Ml, 



fuM/T ■ 



f{(?TflT 



(10) 



which is a perturbative (small magnetic coupling g) loop: 
note that it agrees with the former one due to Dirac con- 
dition e ~ 1/g.) 

If one uses screening masses to get an idea about den- 
sity of electric and magnetic objects, one finds that the 
point at which electric and magnetic masses are equal 
should be close to the E-M equilibrium point we empha- 
sized above. This argument places the equilibrium tem- 
perature somewhere in the region of 



Te=m « (1.2 - 1.5)Tc = 250 - 300 MeV 



(11) 



mJT 




but rather the opposite strongly coupled (correlated) 

screening^. 

Dyons. A very special sector of MQPs arc particles 
with both charges. Because they produce parallel electric 
and magnetic fields, they have nonzero (EB) and thus 
the topological charge. In fact, as shown by Kraan et 
al [25], finite-T instantons can be viewed as being made 
of A^c self-dual dyons: for a very nice AdS/CFT "brane- 
based" construction leading to the same conclusion, see 
[26]. 

Topology is in turn associated with the Dirac zero 
eigenvalues for fcrmions, which can be located and 
counted on the lattice quite accurately. Furthermore, a 
"visualization" of dyons inside lattice gauge field config- 
urations (using variable non-trivial holonomy) has been 
developed into a very sensitive tool [27], revealing multi- 
dyon configurations and their dynamics. One can verify 
that they make a rather dilute but highly correlated sys- 
tems: in fact closed chains of up to 6 dyons of alternat- 
ing charges have been seen. The self-dual dyon density 
and other properties, as well as their relation to instan- 
tons and confinement are summarized in recent paper 
[28]. It is enough to mention only that self-dual dyons, 
like instantons, are electrically screened [29,30] and thus 
rapidly disappear into the QGP at T > Tc. Around Tc 
their density can thus be related to the instanton density 

ndyon ^ Ncriinstantons/T ~ 3/m~^ (13) 

and the mass to the instanton action 

MfiyQyi = T ^ Sinstanton 

/N, ^ (3 - 4)T (14) 

Both are of the order of the density (and the mass) of 
the electric (gluon and quark) quasiparticles at 1.5Tc, 
confirming a suggested E-M equilibrium in this region. 



FIG. 3. Temperature dependence of electric and magnetic 
screening masses according to Nakamura et al [23]. The dot- 
ted line is fitted by the assumption, rrig ~ ff^T. For the 
electric mass, the dashed and solid fines represent LOP and 
HTL resummation results, respectively. 



High-T monopoles. The total pressure related to 
magnetic (3d) sector of the theory and especially the spa- 
tial string tension are other observable related to MQPs 
above T^- for a short recent summary see [24]. Two im- 
portant points made by Korthals-Altes are: (i) MQPs 
must be in the adjoint color representation, to explain 
data on k-strings and magnetic pressure; (ii) there seems 
to be a nontrivial small "diluteness" parameter of the 
MQPs ensemble 



Ml 



{N - l)nM _ l_ 
Ml ~ 20 



(12) 



The fact that screening takes place at distances smaller 
than the average inter-MQP ones is a clear indication 
that screening is not a Debye-type weak coupling one. 



E. Higgs phenomenon in QGP? 

In this subsection we would like to comment, in a 
brief form, on a number of questions which are invari- 
ably asked in connection with Higgs phenomenon and 
monopoles at T > Tc- 

Naively, there is no simple and direct way to apply 
the lessons from supersymmetric theories such as Af=2 
Seiberg-Witten theory to QCD-like setting. The former 



'^If a reader may have doubts that a correlated screening may 
produce such a result, here is an example from the physics 
of the QCD instantons. The typical intcr-instanton distance 
n~^^'^ ~ 1 fm is 5 times larger than the screening length of 
the topological charge Rtop = 1/M{r]') = .2 fm: the corre- 
sponding ratio for monopoles seem to be around 6~^^^ ~ 3. 
In both cases we don't know how exactly the opposite charges 
are correlated: pairs or chains are two obvious possibilities. 
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has scalar fields and flat "moduli space" of possible vacua, 
while the latter has neither scalars nor supersymmetry to 
keep the moduli space flat. 

At finite T the role of Higgs field is delegated to tem- 
poral component Aq of the gauge field: and in fact in glu- 
odynamics there is a spontaneous breaking of the Z{Nc) 
symmetry at T > Tc because the corresponding effec- 
tive action Seff{< Aq >) has Nc discrete degenerate 
minima® . 

Furthermore, the corresponding effective action gets 
small near and large fluctuations in "Higgs VEV" 
< Aq > are seen in lattice configurations; so one may 
think first about a generic case in which it is some (color 
matrix valued) constant in each configuration, to be av- 
eraged with appropriate weight exp[~Seff{< Aq >)] 
later. Thus one may think about an explicit adjoint Higgs 
breaking of the color group, parameterized by Nc — 1 real 
VEVs (e.g. for SU{3) Tr < AqX" >with Gell-Mann di- 
agonal matrices a=3,8). Such breaking makes all giuons 
massive, except the remaining unbroken Nc — 1 U(l)'s 
which remain massless. These remaining U(l)'s are the 
Abelian gauge fields which define magnetic charges of the 
monopoles and their long-range interactions (and electric 
ones, in the case of dyons). 

Finally, the last comment about one lesson from SUSY 
theories which we don't think can be transferred into 
the QCD world: these are the enforced properties of 
monopoles (and many other topological objects like 
branes) which happen to be "BPS states" with their 
Coulomb interactions being exactly cancelled by mass- 
less scalar exchanges. As a result, such objects can often 
"levitate" in SUSY settings. In QCD we however do not 
see or need massless scalars, leaving the usual Coulomb 
and Lorentz forces dominant at large distances. 

II. FEW-BODY PROBLEMS WITH MAGNETIC 
CHARGES 

The simplest few-body system with magnetic charge is 
made of two particles: one has electric charge and the 
other has magnetic charge. In a more general sense we 
should consider them as two dyons, both with nonzero 
electric as well as magnetic charges. This problem has 
been very well studied for many years in both classical 
and quantum mechanics, and it has fascinated the physi- 
cists with many unusual features. See for example [31] 
[32] [33]. 

In such problems one has both electric and magnetic 
fields. We have the electric field from an E-charge (at 
space point Vc) to be 

m = (15) 

\r - Tej'^ 



Fermions will lift this degeneracy, as is well known. 



and magnetic field from a M-charge (at space point rg) 
to be 

^{^=9^^ (16) 

The interaction between one moving dyon {ei,gi) and 
the other (62,172) is given by the Coulomb and 0(v/c) 
Lorentz forces 

-> f 
F12 = (62 • 61 + g2 ■ 91)^ 

+ (62-3i-ff2-ei)^ X 4 + 0(i;Vc') (17) 

with r = r2 — fi. Here we have used the Gaussian units 
in which E and B have the same unit and so are the 
charges e and g. 

As early as in 1904, J. J. Thomson found the even two 
non-moving charges have a nonzero angular momentum 
J carried by rotating electromagnetic field . Indeed, for 
an E-charge and a M-charge (separated by r) as sources, 
it is 

f ^ E X B eg , , 

Jfieid = / d-'xx X = — r (18) 

This angular momentum depends only on the value of 
charges, independent on how far or close they may be. 
Its direction is radial, pointing from the E-charge to the 
M-charge. 



y 




0.0 0.25 0.5 0.75 1.0 



X 

FIG. 4. The trajectory of a dyon in a field of static charge. 

Even earlier, in 1896 Poincare observed that a dyon 
moves in a charge Coulomb field on the surface of a cone, 
as shown in Fig. 4. Their relative motion (angular ro- 
tation and radial bouncing ) is always confined inside a 
cone simply because of the conservation of total angular 
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momentum including the relative rotation and the field's 
angular momentum (18) as well. When getting closer to 
each other the two particles are forced to rotating faster 
thus experiencing effective repelling which makes them 
bouncing radially. Another way to explain the conical 
motion is to notice that a magnetic charge is making 
Larmor circles around the electric field; it shrinks near 
the charge because the fields gets stronger there. 

Quantum mechanics of such two-dyon system has been 
worked out in many details since 1970s, especially the 
many bound states are calculated, see [32,33] for review. 



A. Static Electric Dipole and a Dynamical Monopole 

A very interesting and important few-body problem is 
a magnetic monopole moving in the field of a static elec- 
tric dipole. This is a starting point for studying a " color" - 
electric dipole(quark-anti-quark) surrounded by a gas of 
weakly interacting monopoles which, as we argue, may be 
very much relevant for understanding confinement. Also 
as far as we know, this system seems never been studied 
before. 




FIG. 5. Trajectory of monopole motion in a static electric 
dipole field (with charges at ±1 z) as (left panel) projected 
on x-y plane and (right panel)projected on R-z plane 




FIG. 6. Trajectory of monopole motion in a static electric 
dipole field (with charges at ±1 z) as (left panel) projected 
on x-y plane and (right panel)projected on R-z plane 



Considering quark-anti-quark pair surrounded by a gas 
of weakly interacting monopoles as the scenario for con- 
fining the flux tube just around Tc, one realize that to find 
possible bound states (namely, states with the monopole 
attached around the electric dipole permanently or at 
least for long time before "decaying" away) of such a 
system is potentially a key to understand the source of 
large entropy associated with static quark-anti-quark as 
indicated by lattice. We will discuss this issue in both 
classical and quantum mechanics. 

In classical mechanics, we have the EoM for the 
monopole (of mass m and magnetic charge g) to be the 
following: 



gEd X 



df 
'di 



(19) 



with Ed the electrostatic field from the dipole with ±e 
charges sitting at ±a on z axis: 



E, = e 



az[ 



(20) 



It will be much convenient to work in the cylinder coor- 
dinate (p,(j),z). By running this EoM numerically with 
various initial conditions, we can directly obtain real time 
trajectories of the monopole. 
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FIG. 7. Trajectory of monopole motion in a static electric 
dipole field (with charges at ±1S) as (left paneI)projected 
on x-y plane and (right panel)projected on R-z plane 

A lot of very complicated and very different motions 
have been found, sensitively depending on the initial con- 
ditions. Roughly one may divide these trajectories into 
two categories: "trapping" cases (see Fig. 5 and Fig. 6) 
and "escaping" cases (see Fig. 7). By "trapping" cases we 
mean the monopole starts with ^ a and after a rela- 
tively long time it still remains within distance ~ a from 
the dipole, while in "escaping" cases the monopole begins 
moving further and further away from the dipole after a 
somewhat short time. Due to limited space we show be- 
low only few pictures for both cases. Let's just emphasize 
one particular feature as clearly revealed in Fig. 5: the 
monopole is bouncing back and forth between the two 
electric charges, because of the effective repulsion when 



7 



it is getting close to the charges (as has been explained 
in the charge- monopole motion). We have found many 
such cases which look like two standing charges playing 
E-M "ping-pong" with the monopole. So there are many 
classical bound states for such system, and in principle 
one can scan through the phase space of monopole's ini- 
tial position and momentum to estimate the "trapping" 
states' phase space volume. 

This phenomenon as shown here is dual to the famous 
"magnetic bottle" , a device invented for containment of 
hot electromagnetic plasmas, provided magnetic coils at 
its ends are substituting the electric dipole and a moving 
monopole replaced by the electric charge. 

Now let's turn to the quantum mechanics of such sys- 
tem. One can write down the following Hamiltonian for 
the monopole: 



n 



[p + gAe 
2m 



(21) 



Here A^. is the electric vector potential of the dipole elec- 
tric field, which can be thought of as a dual to the nor- 
mal magnetic vector potential of a magnetic dipole made 
of monopole-anti-monopole. By symmetry argument we 
can require the vector potential as A^ = Af{p,z)4) and 
the monopole wavefunction as ^' = '0(P: z)e^-^'^ with / 
the z-angular-momentum quantum number. Then the 
stationary Schroedinger equation is simplified to be 



-Pz 



2m 



2m p/a h e 



f)? 



(22) 
(23) 



To go further one has to specify a gauge (which is equiv- 
alent to choosing some particular dual Dirac strings for 
the charges) so as to explicitly write down Af. We use 
the gauge which corresponds to the situation with one 
Dirac string going from the positive charge along posi- 
tive z axis to -|-oo and the other going from the negative 
charge along negative z axis to — oo. This gives us: 



z + a 



yjp^ + (z - aY + (2; + a)2 



(24) 



To give an idea of the effective potential we show Fig. 8 
where Veff{p,z,f) = Vejs{\lx^ +y^,z,f) is plotted for 
the x-y plane with z = and / = 0. From the plot we 
can see that there must also be quantum states with the 
monopole bounded within the potential well around the 
dipole for a long time before eventually decaying away. 
Namely one can find states with E = huj -I- iT with F < < 
huj and count how many there are. This issue will be 
separately addressed in details elsewhere. 




FIG. 8. (color online) Quantum mechanical effective po- 
tential for a monopole in a static electric dipole field. See 
text for details. 



B. Flux tubes in a classical gas 

The problem discussed in the previous subsection has 
of course not only "trapped" monopole states emphasized 
above, but also a lot of scattering states for monopoles 
with positive energy. In a simple classical setting it is 
quite clear that in an ensemble (gas) of monopoles, par- 
ticles would scatter off the effective positive potential 
depicted in Fig. 8 and provide effective pressure on it. 
The net result should be nothing but stabilization of the 
electric flux tubes in a gas of monopoles, at a size equi- 
librating gas pressure with the field stress tensor. 

The situation dual to the problem discussed is 
magnetic flux tubes in a gas of electric charges. Since 
this does not need magnetic charges, it takes place in 
electrodynamic plasmas. As a well known example, sta- 
bilized magnetic flux tubes are found in the solar plasma: 
there they extend for millions of kilometers and are seen 
in a modern telescope as a substructure of "solar black 
spots" . 

Although we would discuss classical flux tubes in a sep- 
arate paper, we would like to point out now their role in 
our general picture. We believe those occur at two ex- 
tremes: (i) at T just above Tc, as classical electric flux 
tubes are expelled by dominant MQPs, and (ii) at very 
high T when there are magnetic flux tubes expelled by 
dominant EQPs. In the former case electric flux tubes 
are a natural continuation of quantum electric flux tubes 
which are created by Bose-condensed MQPs in the con- 
fining phase at T < Tc- What we would like to point 
out here is that MQPs expel the electric field lines out 
of the volume they occupy irrespective whether they are 
Bose-condensed or not. 



8 



III. MOLECULAR DYNAMICS WITHOUT 
PERIODIC BOXES 

Molecular dynamics (MD) provides a straightforward 
way to study the various dynamical properties of a clas- 
sical many-body system. The system we are interested in 
is a plasma containing both electric and magnetic (and 
both positive and negative) charges. So in the normal 
convention used by plasma physics community, this is a 
four-component-plasma(FCP). We however would rather 
name it as 2E-2M-plasma to explicitly show its content. 
More widely speaking we may even include one more type 
of particles, namely dyons with both electric and mag- 
netic charges for individual particles, making a 2E-2M- 
4D-plasma. In this paper we will report our results for 
2E-2M-plasma with three different contents: pure elec- 
tric (which reduces to normal TCP) plasma, plasma with 
about one quarter of particles as magnetic charges, and 
plasma with about half of particles as magnetic charges, 
labelled throughout this paper as MOO, M25, M50 respec- 
tively. Comparison among them is expected to give in- 
dications about the role of magnetic charges, especially 
in the transport properties. The microscopic dynamics is 
classical EM, given by Newton's second law together with 
electric Coulomb force (between two E-charges), mag- 
netic Coulomb force (between two M-charges) and 0{v) 
Lorentz force (between one E-charge and one M-charge). 

The routine MD method (as used in GSZ and most 
MD study of usual plasma) is to put desired number of 
particles in a cubic box and then include as many peri- 
odical image boxes (in all three directions) as allowed by 
computing capacity. The summation over images is very 
much time consuming especially for cases with Coulomb 
type long range forces. Also energy conservation is not 
very well preserved after long-time run because of the 
"kick" on particles leaving one truncation boundary and 
entering periodically on the opposite boundary which will 
gradually heat up the system. 

Here it should be emphasized that we have used an 
alternative approach without any periodic boxes. What 
we have done is to simply give all particles certain initial 
conditions and then let them go. It turns out there are 
two different regimes which we deal with separately: 1) 
"plasma in cup" at medium/weak coupling regime, in 
which case we place a sharply rising large radial potential 
barrier at certain radial distance to hold the particles 
inside this "cup"; 2) " self- holding drop" at very strongly 
coupled regime, which means the particles don't fall apart 
into small pieces but behave like a little raindrop and so 
there is no need for a "cup". In this way we are able 
to perform MD easily with thousand particles and can 
conserve energy for less than percent even after really 
long-time run. We will give more technical details about 
our simulations in the second subsection while present 
basic formulae, units and physical parameters in the first 
subsection . 



A. Formulae, units and physical parameters 

For our 2E-2M-plasma, each particle has either electric 

charge or magnetic charge. The E-charges are assigned 
as CjC with Ci randomly and equally given ±1 (ci =0 
for M-charges) and the M-charges are assigned as g^g 
with gi randomly and equally given ±1 (g, = for E- 
charges) too. For a pair of particles their mutual force 
involves three combinations of their charges: = Cj • ej, 
dij = 9i'9j^ and an important new one Kij = ei-gj — gi-Bj. 
In present study we use the same mass m for both types 
of charges. 

The equation of motion for the ith component particle 
is given by: 



dt^ 



C 



K+l' 01- 



rri + - 



cdt 



(25) 



where fji = fi — fj . The first term on RHS is the well- 
known necessary repulsive core without which all classical 
plasma will collapse sooner or later since no quantum 
effect arises at small distance to prevent positive charges 
falling onto negative partners. We choose K = 9 in our 
MD, which is the same as some previous work [11] [34]. 
There is no particular meaning for K = 9 except that we 
want a large value of n which leads to relatively small 
correction (^ l/K) to potential energy between +e and 
— e at and beyond the equilibrium distance. 

To set the units in our numerical study, we use the 
following scaling for length and time (the unit of mass is 
naturally set by particle mass m) 

f = f/ro with ro = {C/e'^)''^ 

i=t/T with T = (mr^/e^)i (26) 

which leads to the dimensionless equation of motion 



1 



+(^)^i^r,. + (^!^)^§xf,. 

e rf- e c rf- dt 



(27) 



With these setting, we have for example: Length = # 
xro, Time = # xr, Frequency = # x^. Velocity = 

2 

# x^. Energy = 7^ x^, etc. All numbers obtained 
from numerics are subjected to association with proper 
dimensional factors in our units. 

Now we still have two dimensionless physical param- 
eters which controls in the above the magnetic-related 
coupling strength: 

1 . .9 = 7 : this parameter characterizes the relative cou- 
pling strength of magnetic to electric sector. In principle 
there is no limitation for it from classical physics. Since 
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we want to focus on the parametric regime which may be 
relevant to sQGP problem near Tc (where electric sector 
gets strongly coupled while magnetic sector gets weak), 
the parameter g is expected to be small, so we will use 
^ = 0.1 in the MD calculation. Naively suppose one has 
a quantum problem with same g = g/e, then by combi- 
nation with the minimum Dirac condition |^ = §-q — t;, 
^ he hc^ 2^ 

one gets = = 1/(25) ~ 5 which is indeed very 
strongly coupled. 

2. 8 = = ^ lJl°. : this parameter tells us how rela- 

tivistic the particles' motion will typically be. The impor- 
tance of this parameter lies in that it controls the strength 
of Lorentz force (/3 • g) between E-M charges. An impor- 
tant observation here is that compared to the Lorentz 
coupling of a pure electric plasma (which is from elec- 
tric current-current) our Lorentz force has only the first 
power of the small parameter /3 and is thus enhanced be- 
cause of the existence of magnetic charges. Since we are 
doing non-rclativistic molecular dynamics, a small value 
of P should be chosen. In the sQGP the typical speed is 
estimated to be about 1/3 of c, in present calculation we 
however will choose jS = 0.1 which on one hand is not far 
from 1/3 and yet on the other hand limits the relativistic 
corrections to be not more than few percent. 

One more physical parameter we should mention here 
is the so called plasma parameter F defined as the ratio 
of average potential kinetic energy (neglecting the sign) 



(28) 



This definition looks a little different from others [11] [34] 
where usual MD study with periodic boxes defines F = 

1^ with a = (3/47rn)^/"^. They use this conveniently 
because in their approach the density is fixed as desired, 
while in our case there is no boxes any more and we use 
the direct ratio which is essentially meaning the same 
thing. The difference and relation between the two F 
values will be further discussed in section VD. 
The plasma parameter F is important in that: 

1. it distinguishes strongly coupled plasma F >> 1 and 
weakly coupled plasma F << 1; 

2. it roughly distinguishes a gas phase F < 1, a liquid 

phase F ~ 10 and a solid (or solid-like) phase F > 100 ; 

3. different types of plasma with the same value of F 
could be compared in order to reveal the dependence of 
macroscopic properties on plasma contents and micro- 
scopic dynamics, and so we will measure properties as a 
function of F. 



<u> 




< U/N > 


<Eh> 




3kBT/2 



tice spacing a = 1.2ro ^. They are given electric charges 
±1 in an alternating way in all 3 directions. Then for the 
M25 (M50) plasma, we randomly pick out 25% (50%) of 
the particles and re-assign them magnetic charges instead 
of electric charges. Then all the particles are randomly 
given initial velocity (for each of the three component) 
Vi,2,3{t = 0) = V* (RANDOMi^) {RAN DOM # is be- 
tween [0,1]) under condition that the total velocity of the 
whole system is zero. For each type of plasma, changing 
the value of V can eventually lead to different equili- 
brated system after certain time. Roughly the larger V 
is, the smaller (lower) the plasma parameter F (temper- 
ature T) will be. The total running time is 2000Ai with 
At = O.It (which is actually our output time step). In 
general it takes ttherm ~ 20 — 30r to equilibrate the sys- 
tem and we start measurements att = 50r. The iteration 
step and accuracy in EoM subroutine is so chosen that 
the energy could be conserved to less than few percent at 
the end of run. As mentioned before, we have two differ- 
ent regimes which we will discuss the details separately 
in the following. 



1. Plasma in cup 

Numerically we found that for about F < 25, the little 
drop we created couldn't hold itself and after some run- 
ning time it will break into a few much smaller pieces, 
which means the surface tension is not large enough to 
maintain the original "big" drop. To confine the particles 
in a finite volume and make them mix up sufficiently, we 
put a radial potential barrier at some cut-distance Rcut 
to make a container holding such plasma: 



V=[B*{r-Rc 



^ e{r - Rc 



(29) 



By choosing B = 5 and L = 11 we make the edge of our 
"cup" a really steep one, thus keeping as many particles 
inside the " cup" as possible at all time because only very 
energetic particles are able to climb up the edge a lit- 
tle and will soon be reflected back. In our simulations we 
have used Rcut = ll''o- In real time of course the number 
of particles confined within the Rcut is always fluctuat- 
ing, so are other macroscopic quantities like energy etc. 
So this system is like a grand-canonical ensemble. For 
this plasma in cup, all the measurements arc made for 
particles inside the cup only (namely with r < Rcut)- 

By looking at the histogram of total number of par- 
ticles at different time points in Fig. 9, one see that the 
system has very good distribution with well-defined aver- 
age ~ 950 and ~ 30 fluctuation width. Another 



B. Details of MD simulations 



In our MD simulations, 1000 particles are initially 
placed on the sites of a 10 x 10 x 10 cubic lattice with lat- 



This value is very close to a = l.lSro which is calculated 
to be the equilibrium value of iVaCZ-like structure under our 
repulsive core. 
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important check is to see if the system is really homo- 
geneous. In Fig. 10 the radial local density n{r) at five 
different time points (from early to very late time) are 
shown, from which it is clear that the density distribution 
is homogeneous and stable enough. The large fluctuation 
at very small r is understandable because one has much 
less particle numbers nAirr^ for small r. One can also 
see that near our cutting edge {Rcut = H) the particle 
density quickly drops down as we want. These observa- 
tions are true in all of our runs, and the number density 
of our cupped plasma at different F is controlled all at 
n 0.17 with negligible variation. These have shown 
that our simulations for cupped plasma is reliable. 
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FIG. 9. (color online) Histogram of total number of par- 
ticles inside Rcut at 1500 different time points. This is an 
example from M25 plasma at F = 0.99. 



kinetic and potential energy at all time points. Clearly 
both distributions make complete sense and so are other 
macroscopic variables which we skip because of limited 
space. Again these justify our "plasma- in-cup" approach. 
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FIG. 11. (color online) Histogram of total kinetic(blue 
square) and potentiaI(red circle)energy inside Rcut at 1500 
different time points. This is an example from M25 plasma 
at F = 0.99. 

The results to be reported in Section IV. and V. are 
all obtained with this method, which cover the F value 
about 0.3 — 14. We want to focus on this region because 
it is most relevant to the sQGP. 

2. Self-holding drop 
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FIG. 10. (color online) Radial density of particles as a func- 
tion of radius r. The five curves are taken from different time 
points. This is an example from M25 plasma at F = 0.99. 

It is also important to check the fluctuation in energy. 
In Fig. 11 we show a typical histogram of fluctuation in 



For about F > 25 we have found our little drop can, 
amazingly, hold itself despite the possible expansion and 
shrinking with considerable amplitude. By mapping the 
particles' coordinates at the end of run we found the 
particles more or less staying around their original po- 
sitions. This very strongly coupled system behaves more 
like a crystal, especially for F 100. In this regime, we 
have found very good collective modes which are shown 
to manifest themselves in the dynamical correlation func- 
tions in a profound way. These results will be reported in 
section VI. It should be pointed out that the self-holding 
region is reached only for pure electric plasma (our MOO 
plasma). For our M25/M50 plasma, with present method 
the largest F that can be achieved (after "cooling" and 
equilibrating scheme) and maintained in a stable way is 
up to ~ 25. The "cooling" method, namely turning on 
a braking force proportional to particle velocity for some 
time and then turning it off, can bring the M25/M50 sys- 
tem down to some instant F ^ 1000 but then the system 
kinetic energy slowly but steadily keeps increasing with 
potential energy getting more negative, the overall effect 
of which eventually increases F back down to few tens. It 
seems indicating the mixture plasma refuses to become 
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solidified even at classical level because of Lorentz type 
force (different from permanent liquid Helium which is 
due to quantum effect). We will leave this issue for fu- 
ture investigation. 



So already from the EoS we've seen considerable differ- 
ence among the three plasma. Since EoS is important 
for dynamical processes, we proceed to study correlation 
functions and transport coefficients in next section, ex- 
pecting to see more differences. 



IV. EQUATION OF STATE 

Before showing the data, we once again emphasize that 
the goal is to compare three types of plasma (MOO, M25, 
and M50) with different E-charges and M-charge concen- 
tration, and all the comparison will be made by plotting 
certain macroscopic properties as a function of plasma 
parameter T. 
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FIG. 12. Temperature T calculated at difFerent plasma pa- 
rameter r in log-log plot for MOO(triangle), M25(square), and 
M50(diamond) plasma respectively, with the three lines from 
linear fitting (see text). 



V. CORRELATION FUNCTIONS AND 
TRANSPORT COEFFICIENTS 

Study of transport coefficients is very important in 
order to understand the experimental discoveries about 
sQGP, such as the very low viscosity and the diffusion of 
heavy quarks. In this section the transport coefficients 
of our three different plasma will be calculated and com- 
pared in order to see the influence of magnetic charges on 
the transport properties. To do that, we will first mea- 
sure certain correlation functions and then relate them 
to the corresponding transport coefficients through the 
Kubo-type formulae, as is usually done in MD works. 

A. Velocity autocorrelation and diffusion constant 

The first correlation function we will study is the ve- 
locity autocorrelation which is defined as: 

1 ^ 

D{r) = — <J2^d^)-^^iO)> (31) 

i—l 

Here r is the correlation time, Vi denotes the velocity of 
the ith particles and the sum is over all particles. The 
average is over thermal ensemble which is done in nu- 
merical program by average over all time points (with 
the number typically of order ^ 1000). 



The first quantity we want to look at is the 
temperature^" dependence on T which is sort of equa- 
tion of state for plasma." In Fig. 12 the EoS for MOO, 
M25, and M50 are compared in log-log plots. Data for 
all three show a linear relation with similar slopes but 
different intercepts. By simple linear fitting we get the 
following parameterized EoS for them: 



AfOO 
M25 
M50 



T = 0.257 /r "-^27 

T = o.i9i/r"-^°^ 
T = 0.125 /r°-^59 



(30) 




^"By temperature T we actually mean fcsT (with the dimen- 
sion of energy in our units) throughout this paper. 

^^Remember in this classical statistical system the kinetic 
energy per particle is given by Ek ~ 3T/2 and total energy 
per particle is E — (1 — r)*i?fe, so the temperature dependence 
on r also gives all information on energy. 
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FIG. 13. (color online) Velocity autocorrelation functions 
D{t) for (from top down at zero time) MOO(black curve), 
M25(green curve) and M50(red curve) plasma, taken at 
r = 1.01, 0.99, 1.00 respectively. 
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In Fig. 13 wc show typical curves for velocity autocor- 
relation function in MOO, M25, and M50 plasma respec- 
tively. A fast damping behavior at small correlation time 
is observed, followed by small fluctuation from random 
noise at longer correlation time. 

The corresponding transport coefficient, namely diffu- 
sion constant, is calculated by the following Kubo for- 
mula 

f-OO 

D= D{T)dT (32) 
Jo 



constant can become few orders of magnitude smaller 
when one changes from weakly coupled gasous regime 
into strongly coupled liquid regime. This qualitative scal- 
ing in coupling is also found from AdS/CFT calculation 
by Casalderrey-Solana and Teaney in [35]. 

Interestingly if one combines (33) and (30), the depen- 
dence of D on T is then obtained: 

MOO : D = \MT^'^^ 
M25 : D = 1A6T"-^^ 
M50 : D = 1.52T°-^2 (34) 
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FIG. 14. Diffusion constant D calculated at different 

plasma parameter T in log- log plot for MOO(triangle), 
M25 (square), and M50 (diamond) plasma respectively, with 
the three Unes from linear fitting (see text). 

In Fig. 14 we plot LogD as a function of LogT for 
M00,M25 and M50 plasma. Approximate linear relation 
is seen for all three, but with visible difference in slopes 
and intercepts. A linear fit gives the following approxi- 
mate functions -D(r): 

MOO : £) = 0.396 /r°-^^2 
M25 : D = 0.342 /r° ™^ 
M50 : £> = 0.273 /r°-^26 ^33) 

At small F < 1 there are considerable differences of the 
three lines. In the physically interesting region F ~ 1 — 10 
the three plasma have visible but not too much differ- 
ence in diffusion constants. The three lines will cross at 
about F ~ 10 and after that deviation from each other 
again grows quickly. The important feature common 
to all three types of plasma as well as to cQGP model 
in [11] is the power-law dropping of diffusion constant 
with increasing coupling strength. We see the diffusion 



B. Stress tensor autocorrelation and shear viscosity 

It is of particular interest to study the shear viscosity 
of our three plasma, as the low viscosity is one of the most 
important discoveries for sQGP from RHIC experiments. 
For this purpose, one can measure the autocorrelation of 
the off-diagonal elements of stress tensor, namely 

^ 1,2,3 

^(r) = ^ < E (T)^fe (0) > (35) 

l<k 

with the stress tensor off-diagonal elements 
^ 1 

i=l ijij 
N N 

= '^m{vi)i{vi)k + '^m{fi}i{ai)k (36) 

i=l i=l 

In the above equations i, j refer to particles while I, k refer 
to components of three-vectors like separation, velocity 

and force. Vij and Fij arc the separation and force from 
particle i to particle j respectively, while fi,Vi, di are the 
position, velocity, acceloration of particle i. The equa- 
valcncc of the two expressions in the second equation is 
discussed in great details in [36] . The V in the first equa- 
tion is the system volume. In Fig. 15 typical plots of r]{T) 
for three plasma are shown. Again the relaxation of real 
correlation is pretty quick and noises dominate the later 
time. With these correlation functions at hand the Kubo 
formula then leads to the following shear viscosity r]: 

POO 

rj= / v{r)dT (37) 
Jo 

In general shear viscosity is a complicated property of 
many-body systems, the value of which depends on many 
factors in a nontrivial way. Roughly, a system with ei- 
ther very small F (like a gas) or very large F (like a solid) 
will have large viscosity while a system in between (like 
a liquid) will have low viscosity with a minimum usually 
in F = 1 ~ 10 (see for example [11] [37]). A qualita- 
tive explanation is that both the particles in a gas and 
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the phonons in a solid can propagate very far (having 
a large mean-free-path) and transfer momenta between 
well-separated parts, thus producing a large viscosity, 
while in a liquid neither particles nor collective modes 
could go far between subsequent scattering, thus making 
momenta transfer very much localized and leading to a 
low viscosity. 
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FIG. 15. (color online) Stress tensor autocorrelation func- 
tions Tj{T) for (from top down at zero time) MOO(black 
curve), M25(green curve) and M50(red curve) plasma, taken 
at r = 1.01, 0.99, 1.00 respectively. 

Now turning to our plasma with magnetic charges, 
since we have a relatively weakly-coupled magnetic sec- 
tor, one may wonder if the magnetic particles will con- 
tribute more to large-distance momenta transfer and 
hence increase the viscosity significantly. We however ar- 
gue that in the opposite, the Lorentz force induced by the 
existence of magnetic charges will somehow confuse parti- 
cles and collective modes, thus helping keep the viscosity 
to be low. Indeed, as shown in Fig. 16, the viscosity goes 
down as increasing concentration of magnetic charges. 
At small r < 1 (in the gas phase) the three curves are 
getting close to each other, but when F increases into 
the liquid region > 1 there is a considerable decrease of 
viscosity in M25 and even more in M50 plasma. The 
M50 with E-charges and M-charges to be 50%-50%, has 
the values of viscosity about half of the pure electric MOO 
plasma at the same F. So, we conclude that the existence 
of magnetic charges may help us to understand the ex- 
tremely low viscosity of sQGP. A rough parametrization 
of the data gives the following viscosity dependence on F 
in the plotted region: 



MOO 
M25 
M50 



1 0.353 



?7 = 0.002 /F^ '''*-f 0.168 /F" 
77 = 0.013 / F i-^^ + 0.105 / F 
7y = 0.096/F"-^°°-)-0.001 • F^-^^ 



(38) 



In all of them the first term is most dominant at very 
small F while the second term becomes important at rel- 
atively large F. We noticed that for M50 there is already 



positive power term of F, which is in accord with the ex- 
pected qualitative feature. Similar terms will appear in 
two other plasma when we will be able to include in the 
fitting more points from large F. 




FIG. 16. Shear viscosity i] calculated at different 
plasma parameter F for MOO(circle), M25 (square), and 
M50(diamond) plasma respectively. 



C. Electric current autocorrelation and conductivity 

The last transport property we study in this paper is 
the electric current autocorrelation and the electric con- 
ductivity. This analysis is only done for pure electric MOO 
plasma since the comparison among MOO, M25 and M50 
(which already have different E-charge concentrations) 
doesn't make much sense. The electric current autocor- 
relation is given by 



a(r) 



1 



N 



N 



3VT 



< (^e,i/,(T)).(^e,«,(0)) > (39) 



1=1 



with Ci the electric charge of the ith particle. And the 
electric conductivity is obtained from Kubo formula as 



a{T)dT 



(40) 



In Fig. 17 we show the typical cr(r) as a function of r 
for two values of F. For this correlation function we do 
notice that even for F not large, the late time correlation 
is not purely noise but still has small oscillation. This is 
not unreasonable since related collective modes like plas- 
mon may develop even for a gas. After integration it 
turns out in the region F ~ 0.3 — 15 the conductivity is 
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scattered between a = 0.101 — 0.141 without clear ten- 
dency, which may indicate the electric current dissipation 
is not sensitive to F in this region. It is very interesting 
to see what will happen to the color-electric conductiv- 
ity (giving information about color charge transport) in 
a Non-Abelian plasma. 



T 




I i \ i I i \ i \ i \ I \ i \ i \ i I i \ i \ i L 
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FIG. 17. (color online) Electric current autocorrelation 
functions cr(r) for pure electric MOO plasma taken at (from top 
down at zero time) F = 6.33(red curve) and F = 14.53(black 
curve) . 



6. We then get the time scale as the inverse of plasmon 
frequency t, = l/a;^ = {j:^f'^ « 1/(4.2T);12 

7. The entropy density is estimated from Stefan- 
Boltzmann limit as s « 0.8 x ^[2 x 8 (7/8) x 2 x 
2x NcX Nf]T^ « 16T3. 

Now let's discuss the value of F. As already mentioned, 
the F given in our MD is the actual ratio of potential to 
kinetc energy, which is measured during the simulation. 

2 

The usually quoted one, defined as F = aws\kBT} ' could 
be considered as a pre-determined 'superficial Gamma'. 
Unfortunately it is not clear how to estimate the actual 
Gamma F of sQGP while the superficial Gamma F is 
obtainable for sQGP, which is F « 2.6 < a^C 2.6. 
So we should try to figure out the superficial Gamma 
in our MD and map the results accordingly. The two 
are different though, they are monotonously related to 
each other, namely when one is large(small) so is the 
other. Since our MD has been done with nX^ — 0.17, the 
{aws)MD ~ 1-12A which means in our MD f = 0.89/r, 
which after combination with (30) will give us the con- 
vertion formula between the two Gamma's. Similar con- 
vertion relation could also be obtained for cQGP in [11] 
from their Fig. 8 though in their case they use superficial 
Gamma as basic parameter and measure potential energy 
from simulation. 



D. Mapping between MD systems and sQGP 

With the MD-obtained empirical formulae for diffusion 
and viscosity, it is of great interest to see what they pre- 
dict for the parameter region corresponding to the sQGP 
experimentally created at RHIC. To do this mapping, 
one has to identify the corresponding physical values of 
basic units (namely mass, length and time) in the desti- 
nation system and then combine dimensionless numbers 
and relations from MD with proper dimensions. Also the 
plasma parameter F should be determined for the desti- 
nation system such that we pick up the MD-predicted 
values of interesting quantities (say, diffusion constant 
and shear viscosity) at exactly the same F value. 

Following similar estimates as in [11], we summarize 
below the relevant quantities of sQGP around l.bTc- 

1. Quasiparticle (quarks and gluons) mass can be esti- 
mated as m ~ 3.0T; 

2. The typical length scale is simply estimated from 
quantum localization to be ~ 1/m « 1/(3. OT); 

3. The electric coupling strength, after averaging over 
different EQPs(quarks and gluons) with their respective 
Casimir, is roughly < a^C >~ 1; 

4. The particle density is roughly given, under the light 
that lattice results have shown the sQGP pressure and 
entropy to reach about 0.8 of Stefan-Boltzmann limit, by 
n « 0.8(0.122x 2 X 8-^0.091 x2x2xA^cXiV/)r3 « 4.2T3; 

5. This density estimation leads to the Wigner-Seitz ra- 
dius aws = i^y^'' ~ 1/(2.6T) « l.lro; 
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FIG. 18. (color online) Plots of Log[l/r;] v.s. Lop[l/L>] for 
three different plasmas. The shaded region is mapped back 
from experimental values, see text. 

With all the above ingredients we are at place to do 
the mapping for interesting transport coefhcients D and 
T] between our MD systems and the sQGP. The mapping 
is a two-way business: one may map the experimentally 
suggested values back into corresponding MD numbers, 
as is done and shown in Fig. 18; or one can use the MD- 
obtained relations to predict the corresponding relations 



This Tp has subtle difference in time scale used in our MD, 
namely the MD time unit r is related to inverse of plasmon 
frequency by r = Tp x (47rnA'^)^''^ « 1.46rp which should be 
taken into account for mapping. 
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of sQGP after conversion of units, as is shown in Fig. 32 
and discussed in the summary part. Here let's focus on 
MD systems in Fig. 18 where a Log[l/ri] v.s. Log[l/D] 
is plotted: data points for all three plasmas fall on a 
universal unit-slope straight line on the left lower part, 
indicating a small F gas limit with diffusion and viscos- 
ity both proportional to mean free path; all three curves 
soon deviate from gas limit at larger F (strong coupling) 
and become flat in the liquid region; the shaded oval 
is obtained by mapping back the following experimental 
values: rj/s ~ 0.1 — 0.3, 2-kTD « 1 — 5, which is clearly 
not close to gas region but near the liquid region, espe- 
cially the one of the M50 curve. More about comparison 
will be given in the summary part at the end of the pa- 
per. 



VI. COLLECTIVE EXCITATIONS AT VERY 
STRONGLY COUPLED REGIME 

In this section we will report interesting results for col- 
lective excitations found at very strongly coupled regime 
(F greater than a few tens) of the pure electric (MOO) 
plasma. The signals of these excitations are extraordi- 
narily clear when F goes to ~ 100 or larger. We came 
to notice these very good modes not in a straight for- 
ward way. Instead, these modes have revealed themselves 
dramatically in some unusual structures of the dynami- 
cal correlation functions and their fourier spectra, which 
we measured first. Only after thinking about possible 
source of these structures we turned to systematic and 
direct measurements for certain collective modes, which 
are found to coincide with correlation functions' struc- 
tures in a distinct manner. As mentioned before, in this 
regime our plasma is like a " self- holding drop" which has 
very different collective motions from a plasma in peri- 
odic boxes: the latter has the familiar phonon modes 
while the former is really like a raindrop, having vibra- 
tion modes like monopole modes, dipole modes, quadru- 
ple modes, etc corresponding to different components of 
density distribution's spherical harmonics. We will dis- 
cuss these modes respectively in more details in the fol- 
lowing. 



A. Monopole modes 

Let's start with the velocity autocorrelation function 
(31) which is supposed to almost vanish (except a little 
random noise) at large correlation time and give conver- 
gent integral to yield diffusion constant, as seen in previ- 
ous section. However when measured in the very strongly 
coupled regime, this correlation function is found to have 
robust oscillating behavior even for very large time which 
of course couldn't be considered as noise, see Fig. 19. By 



looking at the fourier spectrum of it, one immediately 
sees a large and narrow peak at = 0.35 very clearly 
on top of a very broad shoulder structure, as shown in 
Fig. 20. These behaviors are true for F down to about 50. 




Correlation Time 



FIG. 19. (color online) Velocity autocorrelation function 
taken at F = 116.91. 
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FIG. 20. (color online) Fourier transformed spectrum of 
velocity autocorrelation function taken at F = 116.91. 

Now the question is why there will be such peaks in 
velocity autocorrelation. The answer lies in the monopole 
modes, which can be directly measured through simply 
the time dependence of average particle radial position, 
namely: 

1 ^ 

i=l 

In Fig. 21 one sees really nice oscillations lasting for 
near hundred revolutions. In a small movie showing the 
positions of particles at several subsequent time points 
the "drop" in monopole mode looks like a beating heat. 
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This oscillation amplitude decreases slowly indicating a 
nonzero but small width of this monopole mode. Again 
the fourier spectrum gives important information such 
as the characteristic frequency and width of collective 
mode. In Fig. 22 we can see one major narrow peak at 
uji^ — 0.35 together with a few roughly visible but much 
smaller lumps at a; = 0.22, 0.46, 0.70. oj = 0.70 structure 
may be a secondary harmonics of the major peak, and 
seemingly the 0.22 and 0.46 may also be in the same 
series of harmonics with different ranks. 

The important finding we want to point out is the co- 
incidence of U!i^ here with the peak ajf from velocity au- 
tocorrelation function. This result tells us the monopole 
mode, in a form of radial vibration, has nontrivial in- 
fluence on velocity autocorrelation and consequently on 
particle diffusion. 



B. Quadruple modes 

The study of stress tensor autocorrelation in the 
strongly coupled plasma gives us even more interesting 
correspondence between correlation functions and collec- 
tive modes. In Fig. 23 we plot the stress tensor auto- 
correlation (35) as a function of time, and in Fig. 24 its 
fourier spectrum, in which three clear and narrow peaks 
can be seen at lo^ = 0.20, luI^ = 0.40, and uj!^ = 0.45. 
The 0.40 peak, which is the smallest one, may be a sec- 
ondary harmonics of the remarkable 0.20 peak. At F as 
small as about 25, the 0.20 peak is still alive in this cor- 
relation function. Because of the existence of these, the 
correlation function has significant oscillations with large 
amplitude even for very large correlation time. 
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FIG. 21. (color online) Average monopole moment R{t) 
(see text) as a function of time, taken at F = 116.91. 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 
Frequency 

FIG. 22. (color online) Fourier transformed spectrum of av- 
erage monopole moment R(t)(see text), taken at F = 116.91. 




400 600 
Correlation Tiine 



1000 



FIG. 23. (color online) Stress tensor autocorrelation func- 
tion taken at F = 116.91. 
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FIG. 24. (color online) Fourier transformed spectrum of 
stress tensor autocorrelation function taken at F = 116.91. 
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FIG. 25. (color online) Average off-diagonal quadruple 
moment (323 (i) (see text) as a function of time, taken at 
r = 116.91. 
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similar results for the other two. From the figure we can 
see that at the very beginning there is almost no quadru- 
ple mode but its amplitude grows significantly in a time 
interval — 200 during which the monopole mode decays 
down (see Fig. 21). Then after that they persist for long 
time. This indicates that these off-diagonal quadruple 
modes are very robust and somehow "cheap" to excite 
and the energy initially in the monopole modes is prefer- 
ably transferred into the quadruple modes. 

Now when we plot the fourier spectrum of the off- 
diagonal quadruple modes in Fig. 26, amazingly three 
very clear peaks appear at ujf — 0.20, — 0.40, and 
Lu^ = 0.45, which are exactly the same frequencies found 
in the stress tensor autocorrelation. The relative am- 
plitudes among the three peaks are also similar in two 
cases. This is a profound correspondence which means 
the stress tensor correlation and the related transport 
property, namely viscosity, are especially dominated by 
the off-diagonal quadruple modes of the system. This 
type of connection may be universal and one may find 
certain collective excitations for each dynamical correla- 
tion function. 

Before closing this subsection, let's also mention the di- 
agonal quadruple modes which we also studied. This part 
of quadruples could be probed by the following quantity: 

1 ^ 

= nJ2 [^amh)' - \m\'] , ^ = 1,2,3 

(43) 

It has similar behavior as the off-diagonal modes (to 
save space we skip to show the plots) , with peaks in spec- 
trum at Lof =0.12 and lo^ = 0.25 which presumably are 
different ranks in the same harmonic series. These peaks 
however are not seen in the stress tensor autocorrelation, 
which is understandable since the stress tensor autocor- 
relation we study is actually the correlation of stress ten- 
sor's off-diagonal parts which is related to shear viscosity. 
We think these peaks of diagonal quadruple modes must 
be seen in the diagonal parts of stress tensor correlation 
which is related to the bulk viscosity. 



FIG. 26. (color online) Fourier transformed spectrum of av- 
erage off-diagonal quadruple moment (323 (i) (see text), taken 
at r = 116.91. 

To find the source of these, we directly measured the 
off-diagonal quadruple modes by the following probe: 

1 ^ 

Qik{t) ^ -J2^nm{n{t))k ,Z,fc-i,2,3, l^k 

i=l 

(42) 

We have three independent of them, say Q12, Q23, Qsi- 
In Fig. 25 we show one of them as a function of time with 



C. Plasmon modes 

It is clear that our "drop" won't have dipole (or any 
odd-multiple) excitation according to symmetric setting. 
But there can be another type of dipole excitations, 
namely the electric dipole modes, or in a more common 
notion the plasmon modes. These can be probed by 

1 ^ 

{eR)i{t) ^ -Y,e^{nm = 1,2,3 (44) 
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FIG. 27. (color online) Electric current autocorrelation 
function taken at F = 116.91. 
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FIG. 30. (color online) Fourier transformed spectrum of 
average electric dipole moment eR{t) (see text), taken at 
F = 116.91. 



T — ' — I — ' — I — ' — r 




FIG. 28. (color online) Fourier transformed spectrum of 
electric current autocorrelation function taken at F = 116.91. 




FIG. 29. (color online) Average electric dipole moment 
eR{t) (see text) as a function of time, taken at F = 116.91. 



The corresponding correlation function should be the 
electric current autocorrelation (39). To clearly reveal 
these modes we shift our positive charges and negative 
charges with small displacement in opposite directions at 
the initial time which introduces zero dipole but nonzero 
electric dipole. We then plot the electric current autocor- 
relation (Fig. 27) and its fourier spectrum (Fig. 28), and 
the direct electric dipole (Fig. 29) and its fourier spectrum 
(Fig.30) as well. 

Again one see similar behavior in both and find similar 



peak structure at cu'' 



2.0 in both spectra. These 



peaks are large but broad and have some fluctuation, as 
compared with previous peaks, which is understandable 
as the plasmon modes usually have bigger width (larger 
dissipation) than sound modes. These modes seem to 
be present even if the system has only F ^ 10. If one 
calculate the plasmon frequency using the simple formula 

i^nisl for our drop (with n 



l/o3 



1.18A) we 



UJp = 

then get uip = 2.7. This is not far from the observed 2.0 
and the discrepancy must be there because the size of our 
system is only about 10 times the microscopic scale and 
the positive and negative charges in the middle are not 
entirely screening each other as assumed when deriving 
the formula. 



D. Size scaling of the collective modes 

Since the monopole and quadruple modes should be 
sound modes, it is interesting to see how their frequencies 
scale with the system size. For a large enough system 
one expects the sound modes dispersion to he uj = Cgk — 
^Ygk^, namely the mode frequency itself scales linearly 
in k while the width scales quadratically. What we did 
is to change system size to be 10 x 10 x 10, 8 x 8 x 8, 
6x6x6, and 4x4x4, and then look at the change of 
peaks in those monopole and quadruple modes (picking 
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the major peaks wf'^, Ui and ). As demonstrated in 
Fig. 31 where these frequencies are shown as function of 
211 /L with L the system size, the hnear scaUng is very 
weh observed. We obtain the following linear fitting for 
the three modes: 

Monopole: oj^' = 0.610 • k 
Quadruple D : =^ 0.404 • k 
Quadruple N - D: uj'^ ^ 0.329 • k (45) 

The two quadruple modes have similar slope (with the 
diagonal one a little larger) which means they have close 
propagation velocity, while the monopole modes have a 
slope or propagation velocity larger by a factor about 
1.5 — 2. This is reasonable, just like in usual solid the 
longitudinal sound waves have larger velocity than the 
transverse ones. The lines show remarkable consistency 
with the fact that sound modes with infinitely large wave- 
length should have zero frequency. For the width however 
we didn't unambiguously see a regular dependence on L, 
which indicates our systems are not macroscopic enough 
since the width is more sensitive than frequency itself to 
the dissipation effect related to system size. 
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FIG. 31. Peak frequency as a function of 2n / L (with L the 
system size) for monopole modes(circle), diagonal quadruple 
modes (diamond), and off-diagonal quadruple modes(square) 
respectively, with the three lines from linear fitting (see text). 

The last result to mention is about the dependence of 
plasmon modes on system size. Actually we found for all 
four sizes, the plasmon frequencies roughly stay around 
Lo^ = 2.0 without changing, only with the peak structure 
getting worse. Again this is expected since the usual 

plasmon dispersion displays a plateau at ujp = y^^^ 
for k not large. 



VII. SUMMARY AND DISCUSSION 

In summary, we propose to view the finite-T QCD 
as a competition between electrically charged quasiparti- 
cles (EQPs) and the magnetically charged (MQPs). The 
high-T/high density limit is known to be perturbative 
QGP, which is electric-dominated. This implies that 
EQPs are more numerous, with density ~ N^T^, while 
the density of MPQs is - N^T^ /log^{T / Kqcd)- In this 
case the electric coupling is weaker than the magnetic 
e < g. We think that at some intermediate T ~ 300 MeV 
both sectors' couplings and densities are similar, and be- 
low it r < 300 MeT^ the roles are reversed, with dom- 
inant MQPs and electric coupling being stronger than 
magnetic e > g. One of the important consequences of 
this picture is "postconfinement" phenomena and (elec- 
trically) strongly coupled QGP right above the decon- 
finement phase transition T = (200 - 300) MeV. 

Using these ideas as a motivation, we start a pro- 
gram of studies aimed at understanding the many-body 
aspect of matter composed of both of electrically and 
magnetically charged quasiparticles. (Two-body charge- 
monopole and charge-dyon situations have been covered 
in the literature [32,33].) 

Beginning with a 3 body problem: a static E- 
dipole plus a dynamical monopole, we have found evi- 
dences both classically and quantum mechanically that a 
monopole can be bound to an electric dipole, which will 
later be more thoroughly treated in a dilute monopole 
gas scenario and may eventually lead to an explanation 
of large entropy associated with static quark-anti-quark 
just above Tc as indicated by lattice data. 

We have then used molecular dynamics to do the first 
systematic study of a plasma with both electric and mag- 
netic charges. Two regimes have been separately studied: 

1. In the weak and medium coupling regime (F < 25), 
which is most relevant to sQGP problem, the equa- 
tion of state, the correlation functions and the transport 
coefficients have been evaluated and compared among 
plasma with different magnetic contents. Most interest- 
ingly we found by increasing the concentration of mag- 
netic charges to about 50% we get a factor 2 down for 
viscosity, which is particularly important in view of ex- 
plaining surprisingly low viscosity of sQGP as observed 
at RHIC. 

2. In the very strongly coupled regime (F > 25), very 
interesting correspondence between correlation functions 
and collective excitations have been revealed. The 
monopole modes are found to induce long time oscilla- 
tion in velocity autocorrelation, while the off-diagonal 
quadruple modes are shown to have profound dominance 
on the nontrivial oscillating behavior of stress tensor (off- 
diagonal) autocorrelation. In a similar manner the plas- 
mon modes are connected to electric current autocorrela- 
tion. So with such relation, studying one side also gives 
us information about the other side. The sound modes 
nature of such collective excitations has been dcmon- 
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strated by studying the scaling property of their peak 
frequencies with system size. 

Finally, we would like to compare our results with 
those obtained using the AdS/CFT correspondence and 
also with empirical data about sQGP from RHIC exper- 
iments. Those are summarized in Fig. 32, as a log-log 
plot of properly normalized dimensionless (heavy quark) 
diffusion constant and viscosity. 



AdS/CFT will only become quantitative with time: this 
figure is just the first attempt to get together all three 
ingredients of the broad picture. 
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FIG. 32. (color online) Plots of Lop[l/(7?/s)] v.s. 
Lo(7[l/(27rT_D)] including results from our MD simulations, 
the Ads/CFT calculations, the weakly coupled CFT calcula- 
tions, as compared with experimental values, see text. 

The dashed curve in the left lower corner is for A/'=4 
SUSY YM theory in weak coupling, where viscosity is 
from [41] and diffusion constant from [40]. The curve 
has a slope of one on this plot, as in weak coupling both 
quantities are proportional to the same mean free path. 
As one can easily see, weak coupling results are quite 
far from empirical data from RHIC, shown by a gray 
oval in the right upper corner. Viscosity estimates follow 
from deviations of the elliptic flow at large pt from hydro 
predictions [38], and diffusion constants are estimated 
from Raa and elliptic flow of charm [39] . 

The curve for strong-couphng AdS/CFT results (vis- 
cosity according to [42] with 0(A~'^/^) correction, diffu- 
sion constant from [35]), shown by upper dashed line, 
is on the other hand going right through the empirical 
region. At infinite coupling this curve reaches s/77 = 47r 
which is conjectured to be the ever possible upper bound. 

Our results - three solid lines on the right - corre- 
spond to our calculations with different EQPs/MQPs ra- 
tio. They are close to the empirical region, especially the 
version with the equal mixture of EQPs and MQPs. 

Let us end with a warning, that the empirical data, 
the mapping from classical system to sQGP and the re- 
lation between QCD and the M =4 SUSY YM used in 
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